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Abstract 

Nonequilibrium molecular dynamics simulations is used to study the 
motion of a Cqo molecule on a graphene sheet subjected to a temperature 
gradient. The C(,o molecule is actuated and moves along the system while 
it just randomly dances along the perpendicular direction. Increasing the 
temperature gradient increases the directed velocity of Cqq. It is found 
that the free energy decreases as the Ceo molecule moves toward the cold 
end. The driving mechanism based on the temperature gradient suggests 
the construction of nanoscale graphene-based motors. 



1 Introduction 

Since graphene has been discovered [JJ , many properties of this two-dimensional 
material have been studied both experimentally [5] and theoretically [3J H] . In a 
recent experimental research the dynamic of light atoms deposited on a single- 
layer garphene has been studied by the mean of transmission electron microscopy 
(TEM) technic [5]. The two-dimensional structure of graphene suggests the 
possibility of motion with just two degrees of freedom. The free energy surface 
for a particle moving above a graphene sheet explains different motion-related 
phenomena at nanoscale as well as the various directed motions on the carbon 
nanotube-based motors [51 [7] • 

A net motion can be obtained from a nanoscale system subjected to a ther- 
mal gradient [HI IS]- Recently the motion of an experimentally designed nanoscale 
motor consisting of a capsule-like carbon nanotube inside a host carbon nan- 
otube has been explained successfully with molecular dynamics (MD) simula- 
tions [10]. The capsule travels back and forth between both ends of the host 
carbon nanotube along the axial direction. Barreiro et al. have designed an 
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artificial nanofabricated motor in which one short carbon nanotube travels rel- 
ative to another coaxial carbon nanotube 1.7; . This motion is actuated by a 
thermal gradient as high as 1 K nm _1 applied to the ends of the coaxial carbon 
nanotubes. 

Since graphene has a very high thermal conductivity (3000-5000 W K _1 m _1 [TTJ 
[T2"]). as high as diamond and carbon nanotubes [T3J [TJl [TS1 [TB] , it is a good can- 
didate for heat transferring designs in nano-electromechanical systems. Because 
of strong covalent bonds in graphene, thermal lattice conduction dominates the 
electrons contribution [TT]. Recently Yang et al. [T7] have studied the thermal 
conductivity and thermal rectification of trapezoidal and rectangular graphene 
nanoribbons and found a significant thermal rectification effect in asymmetric 
graphene ribbons [i~7] . 

Here we study the motion of a nanoscale object, e.g. Ceo, on a graphene 
sheet, in the presence of a temperature gradient. We show that the graphene is 
a good two-dimensional substrate for thermal actuation due to its high thermal 
conductivity, however as it is expected, in the absence of a thermal gradient, the 
6*60 molecule randomly diffuses on the graphene sheet [5] . The average velocity 
along the temperature gradient direction and the free energy change throughout 
the system are calculated. 

This paper is organized as follows. In Sec. 2 we will introduce the atomistic 
model and the simulation method. Sec. 3 contains the main results including 
those for the produced temperature gradient, the trajectory of the Ceo molecule 
over the graphene sheet and the free energy change. A brief summary and 
conclusions are included in Sec. 4. 

2 The model and method 

The system was composed of a graphene sheet as a substrate, with dimensions 
L x x L y — 70x5nm 2 and a Ceo molecule above the sheet. The graphene sheet 
with N = 14400 carbon atoms was divided into 12 equal rectangular segments. 
Each segment with Ni = 1200 carbon atoms was arranged in 40 atomic rows 
(along armchair or x-direction) . Each row has 30 atoms which were arranged 
along zigzag direction. The system was equilibrated for 300 ps at T = 300 K 
before the temperature gradient was applied. Once the system was equilibrated, 
the first (hot spot) and last (cold end) segments of the graphene sheet were kept 
at T/j and T c , respectively. A temperature gradient between the two ends was 
then produced, i.e. (T/, — T c )/L x (top panel of Fig. [T|). To make the model 
more efficient and prevent crumpling (see Fig. [2£a)) of the ends, we fixed the 
z-components of the first atomic row in the first segment and the last row of 
the last segment (see Fig. [2^b)). 

We carried out MD simulations employing two types of the interatomic po- 
tentials: 1) the covalent bonds between the carbon atoms in the graphene sheet 
and in the Ceo molecule are described by Brenner potential [T^ and 2) the non- 
bonded Van der Waals interactions between the graphene atoms and those of the 
Cqo molecule. Brenner potential has been parameterized to model sp 2 covalent 
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Figure 1: (Color online) Top: The model which shows applied temperature 
gradient along x-direction. Bottom: Produced four temperature gradients after 
450 ps of a nonequilibrium molecular dynamics simulation. Delta symbols are 
related to AT = 100 K, and for square symbols AT = 80 K, circle symbols 
AT = 60 K and for gradient symbols AT = 40 K, respectively. 
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Figure 2: (Color online) (a) A snapshot of a not-constrained system shows the 
effect of the ends crumpling, (b) Fixing z-components of the atoms at the first 
and last rows (see the text) prevents the crumpling. The black dots show a 
typical trajectory of the motion of a Ceo molecule over the graphene sheet from 
the hot spot toward the cold spot. 
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bonds in the graphene, carbon nanotube and Ceo structures. For the non- 
bonded potential a Lennard- Jones (LJ) potential gives reasonable results [19]. 
Here we choose the LJ parameters as e = 2.413 meV and a = 3.4 A [50] which 
represent the depth and range of the LJ potential energy, respectively. Note 
that the LJ potential is a simple and commonly used potential for modeling the 
interaction between carbon nanostructures [5TJ [55J . The equations of motion 
were integrated using a velocity- Verlet algorithm with a time step At = 0.5 fs. 
The temperature of the hot end (Th) and the cold end (T c ) were held constant 
by a Nose-Hoover thermostat. The temperature of the inner segments were not 
controlled by the thermostat. Periodic boundary condition was applied only in 
the y-direction. Due to the applied temperature gradient, the system can no 
longer be described by equilibrium methods and we shall thus employ nonequi- 
librium molecular dynamics simulations. A temperature gradient were produced 
across the system during 450 ps and a stationary state was established. The Cqq 
molecule was put above the second segment at z cm = 8 A (here index cm refer 
to the center of mass). During the production runs of 750 ps both the x and y 
positions of the center of mass of the Cqq molecule were recorded. Moreover the 
temperature of each segment was calculated by measuring the total kinetic en- 
ergy of that segment. In our simulations a typical value for the relative standard 
deviation of the total energy of the extended system is about 3.5 x 10 -5 . A full 
simulation run takes about 50 h CPU time on a 3.2 GHz Pentium IV processor 
with 4 GB RAM. 

In order to compute the change in the free energy, one can employ the 
commonly used thermodynamic integration and perturbation methods |23j . 
A good estimation for the absolute value of the free energy requires sampling 
the whole phase space which is not feasible. Jarzynski's method removes this 
difficulty for nonequilibrium simulations [23]. There is an equality between the 
change of the free energy F and the work W applied on the system (here the 

C 6 o molecule) [24] 

AF = -ft^ 1 lnexp(— /SW), (1) 

where /3 = l/fc^T (T is the temperature in each segment), and the average 
is taken over different configurations with different initial conditions. In fact, 
Eq. ([T]) connects the change of the free energy (between two equilibrium state) 
and the applied work on the system in a nonequilibrium process. 

3 Results and discussion 

3.1 Producing temperature gradient 

The temperature profiles for different temperature gradients with AT = 40, 60, 80, 100 K 
are shown in Fig. [T] (bottom panel). In this figure, the local temperatures of 
each segment which were obtained by averaging over 500 data are indicated 
by symbols. Corresponding error bars indicate the statistical errors and are in 
the range 4-6 K. Notice that the temperature profiles are nonlinear which is a 
commonly observed behavior in nonequilibrium molecular dynamics simulation 
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Figure 3: (Color online) (a, b) Time series of x and y components of the center 
of mass of a Cgo molecule over a monolayer graphene subjected to a temperature 
gradient. Thick curves show the average curves. The inset shows the variation 
of the velocity versus temperature gradient, (c) The trajectory of the motion of 
a Ceo molecule (in x — y plane) moving over the monolayer graphene, averaged 
over six simulations with different initial conditions. 




Figure 4: (Color online) Mean square displacement for y-component of the 
center of mass of Ceo molecule over the graphene. The inset shows the velocity 
autocorrelation function for the y-componcnt of the motion as a function of 
time. 
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Figure 5: (Color online) (a) Total work performed on Ceo molecules and (b) the 
change of free energy during Ceo motion from the hot end to the cold end. The 
thick curve is the average over six simulations with different initial conditions. 



of thermal conductivity [HJ [T51 [TB] . It is a consequence of the strong phonon 
scattering caused by the heat source or heat sink and can be explained by partly 
diffusive and partly ballistic energy transport along the system |14l [25] . Also, 
ripples in the graphene (Fig. HJb)) are the other important source for phonon 
scattering [26, 27 . Therefore, these mechanisms cause a non-linear temperature 
profile in the middle segments. 

3.2 Trajectory in the x — y plane 

The graphene ribbon is a two-dimensional way for the motion of the Cgo molecule 
along it. Ceo attempts to find its equilibrium state and looks for the local min- 
imum of the free energy. In other words, in the presence of a temperature 
gradient, the phonon waves created in the hot spot travel through the system, 
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interact and transfer momentum to the Cgo molecule which results a net mo- 
tion [7]. The time series for the x and y coordinates, separately, are shown in 
Figs. (2Ha,b). The trajectories of the six molecules (from the six different 
simulations) in x — y plane over a time interval of 750 ps are depicted in Fig.(3Jc). 
In all three cases the thick curve is the average of the six others. 

The motion along the ir-direction is almost a linear uniform motion with 
velocities v x = 53, 56, 59 and 63m/s for AT = 40, 60, 80, 100 K respectively 
(see inset of Fig. (3Ja)). There are two computational methods that can be used 
to show that the motion along the y-direction is diffusive (Fig. [3fc)). First, one 
can use the Einstein relation to find the diffusion constant D y by measuring 
the slope of the mean square displacement (MSD) of the y-components of the 
position, i.e. ((y(t) — y(0)) 2 ). To show that the motion along the y-direction 
is driftless, we look at the MSD of y-componcnt of the position of the Cqo 
molecule. For a driftless motion the MSD should grows linearly with time, i.e. < 
(y(t)-y(O)) 2 >= 2D y t. We depict the MSD in Fig. g]for a long time simulation. 
As we see from this figure the slope of the MSD is almost one (after equilibration) 
which is a signature of diffusive regime. Alternatively one can use the Green- 
Kubo relation which makes use of the velocity autocovariance function ( VACF) . 
More specifically, one can take the integral over the autocovariance of the y- 
components of the position of the velocity of Ceo as D y — J (v y {0)v y {T))d,T, 
which gives the diffusion constant of random walk motion along the y-axis. The 
inset of Fig. 2] shows the VACF versus time. We have tested both methods 
in order to compute D y . The result for the diffusion constant is about 4 x 
10- 9 m 2 /s. 

The directed motion resulted from the temperature gradient along the graphene 
sheet, provides a nanoscale motor for the material transferring. Notice that 
in the absence of temperature gradient, when the system is equilibrated at T 
= 300 K, we found a diffusive motion on the graphene sheet [6]. 

3.3 Free energy reduction 

Figures [5ja,b) show the variation of the total work and the changes in the free 
energy with time, respectively. The Cgo loses on average almost 12 eV of free 
energy after 750 ps during the motion towards the cooler region. Note that here, 
the condition W > AF is always satisfied which is a well known criterion for a 
nonequilibrium (irreversible) thermodynamics evolution |24j . This method for 
calculating the change of the free energy is a well known method in computa- 
tional soft condensed matter [28] but it turns out to be useful also for studying 
various physical properties of graphene, particularly investigating the stability 
of new designs for nanoscale molecular devices as studied here. 

4 Conclusions 

In summary, Cgo moves directly toward the cold end of the graphene sheet when 
a temperature gradient is applied along armchair direction. It is found that the 
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lateral motion (along zigzag direction) is a diffusive motion. The reduction of 
the free energy of the system along the molecule motion is an indicative of the 
drifted motion. Comparing the free energy difference with the work performed 
on Ceo shows that the process is indeed thermodynamically irreversible. The 
proposed mechanism for driving nanoparticles on a graphene sheet may be used 
in the design of novel nanoscale motors. 
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